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NEUTRON  FLUX  DISTRIBUTIONS  IN  MULTIPLE  REGION  REACTORS 


-  ABSTRACT  - 

In  this  report  the  problem  is  considered  of  determining  neutron  flux 
distributions  in  reactors  consisting  of  discrete  regions  of  substantially  dif¬ 
ferent  nuclear  properties o 


Starting  with  these  equa¬ 
tions  the  problem  is  considered  of  determining  the  flux  in  one  of  the  energy 
groups  given  the  spatial  flux  distribution  of  the  next  higher  energy  group0 
The  problem  is  handled  by  converting  the  differential  equation  with  source  term 
to  an  integral  equation^  under  appropriate  conditions  the  integral  equation 
may  be  solved  by  a  numerical  iterative  procedure «  The  conditions  under  which 
this  may  be  done  are  discussed  and  formulas  are  developed  for  placing  bounds 
on  the  errors o  In  addition  the  technique  of  obtaining  flux  distributions  in 
complicated  configurations  by  additive  and  multiplicative  superposition  of 
flux  distributions  of  simpler  configurations  is  considered*,  Expressions  are 
derived  for  estimating  the  errors  made  in  using  these  methods o  A  number  of 
simple  examples  are  presentedo 
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SYMBOLS 


A,  B 


D. 

i 


"f 

G 

Ri 

S. 

i 


1 

r« 


see  page  16 

neutron  diffusion  constant  in  i*th  medium 
thermal  neutron  group  diffusion  constant 
fast  neutron  group  diffusion  constant 
Greenes  function 
the  i8  th  region 

the  surface  of  the  iBth  region 
material  multiplication  constant 
resonance  escape  probability 

source  strength  of  neutrons  in  the  i!th  region 
position  vectors 


a,  p,  Y 

€n 

Ki 

2 

9 

^s 

<Pf 

V 

d/dn 


see  page  18 

error  in  the  n*th  iterate 

(Xj/D.)1/2 

neutron  absorption  cross  section 

thermal  neutron  group  absorption  cross  section 

fast  neutron  group  absorption  cross  section 

neutron  flux 

thermal  neutron  flux 

fast  neutron  flux 

gradient  operator 

derivative  in  the  direction  of  the  normal 
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NEUTRON  FLUX  DISTRIBUTIONS  IN  MULTIPLE  REGION  REACTORS 


Sylvan  Nallach 


INTRODUCTION 


The  studies  reported  in  this  paper  were  carried  out  to  develop  methods 
for  determining  the  flux  distributions  in  reactors  made  up  of  discrete  regions  of 
substantially  different  nuclear  properties o 

A  procedure  for 

calculating  the  flux  distribution  by  an  iterative  technique  is  presented  here  as 
well  as  certain  analytical  results  pertaining  to  the  general  problem0  This  work 
was  carried  out  at  the  Westinghouse  Atomic  Power  Division  under  Atomic  Energy 
Commission  Contract  AT-ll-l-GEN-14 „ 
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The  problem  of  determining  the  thermal  neutron  flux  for  a  reactor 
having  discrete  regions  has  been  considered  within  the  framework  of  multi-group 
diffusion  theory,,  A  two-group  model  is  used„ 

A  very  considerable  simplification  of  the  problem  can  be  made  by  assuming  in 
the  two-group  case  that  the  fast  flux  spatial  distribution  is  known  or  may  be 
satisfactorily  approximated <>  For  the  multi-group  case  a  similar  remark  applies 
to  the  flux  in  the  group  next  higher  in  energy  than  the  one  under  consideration,, 
From  a  mathematical  point  of  view  the  proolem  then  becomes  one  of  finding  ap¬ 
propriate  solutions  to  the  Helmholtz  equation  with  sources  in  multiple  regions 
with  boundary  conditions  given  for  the  function  and  possibly  its  gradient,, 

The  problem  outlined  may  be  attacked  by  direct  numerical  solution  of 
the  differential  equation  with  boundary  conditions  „  This  is  the  technique 
-^Viich  has  been  used  successfully  by  Garabedian  and  Schiffo  Tney  have  applied 
the  relaxation  method  of  Southwell  to  determine  the  flux  distributions  at  and 
about  cross  shaped  water  channels „  Another  approach,  the  one  developed  in  this 
report,  is  to  convert  the  differential  equation  and  boundary  conditions  to  an 
integral  equation,.  The  integral  equation  lends  itself  well  both  to  analysis 
and  to  solution  by  numerical  procedures „  No  comparison  of  the  two  techniques, 
insofar  as  numerical  solution  is  involved,  has  been  attempted,.  It  is  likely 
that  both  methods  have  their  respective  advantages  which  may  be  significant  in 
specific  situations o  The  integral  equation  approach  is,  as  will  be  seen. 


A  Ho  Lc,  Garabedian  and  R„  R„  Schiff, 

WAPD-FM-89. 
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amenable  to  analysis  and  permits  bounds  to  be  set  on  the  errors  involved  in  the 
numerical  solution*.  Further  it  appears  that  fairly  good  results  can  be  obtained 
with  it  starting  with  poor  initial  approximations 

Part  I  of  this  report  is  concerned  with  the  derivation  of  the  integral 
equation  and  a  discussion  of  its  form*,  An  iterative  procedure  for  solution  of 
the  equation  is  developed  and  the  conditions  under  which  the  iterative  proce¬ 
dure  converges  is  studied 0  It  is  shown 

that  the  iterative  procedure  converges  for  a  convex  water  hole0 
Bounds  for  the  errors  in  the  iterates  are  developed  under  these  convergent  con¬ 
ditions  o  It  is  also  shown  that  the  iteration  will  fail  to  converge  -under  cer¬ 
tain  circumstances o  Simple  examples  of  a  slab  water  hole  and  a  cylindrical 
water  hole  in  an  infinite  reactor  are  solved  to  illustrate  the  technique <> 

Finally  a  cross  water  hole  is  considered  to  illustrate  the  procedure  for  a  non- 
convex  region*. 

Part  II  consists  of  a  discussion  of  the  technique  of  additive  super¬ 
position  to  obtain  flux  distributions  by  addition  of  the  distributions  for 
simple  regions  and  the  technique  of  multiplicative  superposition  to  obtain  flux 
distributions  for  finite  reactors o  As  an  example  of  additive  superposition  a 
cross  water  channel  is  solved  by  superposing  cylindrical  water  holes 0  Expres¬ 
sions  for  the  errors  introduced  in  superposition  are  developed  for  both  additive 
and  multiplicative  superposition 0  Several  additional  examples  in  slab  geometry 
are  included  to  illustrate  the  methods o 
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I -  THE  INTEGRAL  EQUATION  FOR  DIFFUSING  MEDIA 

This  part,  of  the  report  is  devoted  to  the  derivation  of  the  integral 
equation,  a  description  of  its  use,  and  its  application  to  some  simple  problems . 
It  will  be  shown  how  bounds  for  the  flux  can  be  estimated  and  sufficient  condi¬ 
tions  for  the  convergence  of  successive  approximations  will  be  derivedo 

1»  The  Two  Group  Diffusion  Equations 

The  steady  state  two  group  diffusion  equations  are 

(I  D  V  DgV  cps  -  2S(PS  +  P  If<Pf  =  0 

(1.2)  VDfVcpf~  Zf(pf+|2s(ps  =  0  . 

Suppose  we  are  interested  not  in  criticality,  but  only  in  the  thermal  neutron 
flux  distribution0  Then,  rather  than  attempt  to  solve  a  complicated  charac¬ 
teristic  value  problem,  it  is  frequently  sufficient  to  guess  at  the  form  of  <pf 
and  to  solve  (Id)  for  <p  „  The  slow  flux  thus  obtained  is  often  a  very  good 
approximation  to  the  correct  flux0*  Instead  of  (Iol)  and  (Io2),  let  us 
therefore  consider  for  the  present 

(I03)  VDVcp  -  £<p  +  q  =  0  , 

where  we  are  now  interested  only  in  the  thermal  neutron  flux  and  q  =  p  Yf(?f  is 
the  source  term  for  thermal  neutrons  o  01  course  cp  can  be  any  neutron  energy 
group  of  interest  and  q  the  appropriate  source  tenn0  It  is  equation  (Io3) 


A  If  desired,  the  assumed  values  for  (p„  can  be  checked  by  putting  the  calculated 
(ps  in  equation  (l02)  and  solving  for  (p^o  If  the  iteration  is  continued,  and 
if  the  solutions  obtained  are  suitably  normalized,  the  iterations  can  be  ex¬ 
pected  to  converge  to  the  fundamental  modes  of  the  fast  and  thermal  fluxes <> 
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to  which  Garabedian  and  Schiff  have  applied  the  relaxation  technique  and  which 
we  shall  now  convert  to  an  integral  equation 0  Frequently  equation  (lo3)  is 
referred  to  as  a  one  group  theory  although  it  is  much  more  closely  related  to 
two  group  diffusion  theory  than  to  the  one  group  theory  from  which  is  derived 
the  equation  VhVcp  +  (k—l)  2(p  —  Oo  It  would  be  apt  to  say  that  (lo3)  is 
based  on  a  modified  two  group  theory,, 

2  o  The  Modified  Two  Group  Integral  Equation 

In  this  section  the  integral  equation  defining  the  thermal  flux  dis¬ 
tribution  in  a  region  consisting  of  three  homogeneous  diffusing  media  will  be 
derived  from  equation  (lo5)o  The  usual  boundary  conditions,  i^e,  continuity  of 
flux  and  of  current  at  interfaces,  will  apply □  On  finite  boundaries  cp  =  0;  if 
the  medium  extends  to  infinity  it  is  required  that  cp  remain  bounded 0  The  latter 

/QG 

cp  >  where  G  is  a 

Greenes  function  to  be  defined  beloWo  Consider  some  arrangement  of  the  media 
such  as  is  shown  in  Figure  1. 

Let  us  rewrite  equation  (lo3) 

(I„4)  VDVcp  -  ^cp  =  (£  -  ^1)(p  -  q  , 


ana  consider  also 

(1,5)  DjV2G  -2^  =  0  • 

The  solution  G  of  (X«5)  is  the  Greenes  function  appropriate  to  the  total 
volume  under  consideration*  and  the  boundary  conditions <>  It  is  the  flux  dis¬ 
tribution  created  by  a  point  source  of  neutrons o  If  the  region  is  finite,  then 
G  vanishes  on  the  outermost  boundary,  S;  if  the  region  extends  to  infinity, 
then  G  tends  to  zero  at  infinity,.  More  precisely,  in  the  two  dimensional  case, 
G(r,riI)  is  the  solution  of  (lo5)  which  vanishes  on  the  boundary,  S,  which  has  a 
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logarithmic  singularity  at  r9  and  -which  is  normalized  so  that 

f  o  = 

If  the  medium  extends  to  infinity  in  all  directions  then  G  =  Kq(*i|  r~rff  |  )/2tt* 
where  is  the  zeroth  order  Bessel8 s  function  of  the  second  kind  with  imaginary 
argument* 

Let  us  call  region  the  external  region  and  regions  and  Rg  in¬ 
ternal  regions*  Let  the  point  of  interest  be  an  internal  point,  say  at  r  in 
R^o  Multiply  equation  (1*4)  by  G,  equation  (1*5)  by  cp,  subtract,  and  integrate 
over  all  three  regions  with  the  exception  of  a  sphere,  2,  of  radius  e  around 
the  point  r*  Transform  the  integrals  by  Greenes  theorem  and  let  €  —  0* 

Observe  that 


e—  0  S 


lim 
€  — ►  0 


There  results,  in  view  of  the  boundary  conditions  and  the  continuity  of  G 
and  its  gradient, 

(1.6)  -D,(p  +  (D,-Dp)  /  V  (p-VG  +  (D,-D-)  lim  J  V( p»VG 

=  /  o/(2  -I,)<p  -  i7  • 

R  ■L 

1+2+3 


/  VrVo  =  /  >pg-  / 

sz  h 


/  fS-  / 


f  Zi 

L ”T‘ 
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It  is  evident  from  the  symmetry  in  the  right  hand  side  of  (1.7)  that 
exactly  the  same  expression  will  be  obtained  for  D^cp  if  the  point  r  is  in  region 
R  o  It  is  also  clear  how  additional  regions  can  be  included.  Perhaps  not  so 
obvious  is  the  fact  that  for  points  in  region  R^,  D^cp  also  is  given  by  the  right 
side  of  (1.7)0  It  must  be  remembered*  however*  that  the  singularity  of  G  is  at 
the  point  r  of  interest*  and  that  the  contributions  of  the  surface  integrals 
undergo  discontinuities  as  the  point  r  crosses  an  interface  between  media 0  Thus 
the  content,  if  not  the  form,  of  the  right  hand  side  of  (1.7)  is  different  in  the 
three  caseso 

Equation  (1.7)  is  an  inhomogeneous  linear  integral  equation  for  the 
flux  cp.  It  should  be  noted  that  the  integrals  containing  cp  involve  only  the 
internal  regions  Rg  and  R^.  Thus  the  integral  equation  can  be  solved  for  the 
flux  in  regions  Rg  and  R^  without  any  consideration  of  the  flux  in  R^,>  After 
the  flux  in  Rg  and  R^  has  been  determined,  the  flux  in  region  R^  can  be  ob¬ 
tained  by  integration o  This  fact  is  important  in  solving  (lo7)  by  iteration. 


Ik 


It  is  necessary  only  to  assume  trial  values-  for  the  flux  and  to  iterate  in 
regions  and  R„  in  order  to  determine  the  flux  there <, 

In  solving  (Xo7)  by  iteration  the  flux  on  the  boundaries  of  regions 

R  and  R  can  be  found  in  two  ways,  according  as  the  boundary  points  are  con- 
il  o 

sidered  internal  points  or  external  points o  In  some  cases  the  flux  at  external 
points,  obtained  by  one  iteration,  is  very  much  less  sensitive  to  the  starting 
values  than  is  the  internal  fluxQ  Hence  it  may  sometimes  be  advantageous  to 
estimate  starting  values  for  the  iteration  in  regions  Rr>  R^j  iterate  once  to 
obtain  improved  values  on  the  boundaries  of  R^  and  R^,  and  then  to  make  new  es¬ 
timates  of  the  flux  in  Rg  and  R3  from  these  boundary  values o 

The  integral  equation  (I07)  divides  the  contributions  to  (p  neatly 
into  three  parts— the  first  part  caused  by  the  differences  in  the  diffusion 
constants,  the  second  part  caused  by  the  differences  in  the  /cD s,  and  the  third 
part  caused  by  the  differences  in  the  source  terms o  This  fact  can  be  very 
useful o  In  a  given  problem  it  will  be  clear  how  to  alter  the  values  of  the  con¬ 
stants  so  as  to  obtain  simpler  integral  equations  whose  solutions  majorize  and 
minor ize  the  flux  sought 0  It  is  particularly  advantageous,  from  the  point  of 
view  of  computation,  to  eliminate  the  surface  integrals a  The  extreme  example 
of  this  device  is  the  majorant  provided  by  the  maximum  value  of  q^/  £  ^  and  the 
minorant  provided  by  the  minimum  value  of  iv.ec  the  largest  and  smallest 

asymptotic  values  of  the  fluxes 0 

3 o  Formal  Solution  by  Iteration 

In  this  section  the  solution  of  the  ^ntegral  equation  (lo7)  by  the 
method  of  successive  approximations  will  be  discussed  in  a  purely  formal  mannero 
Questions  of  convergence  are  deferred  to  sections  4  and  7o  We  consider  only 
the  physically  interesting  case  in  which  D . ,  2.,  q.  are  non-negative  and  all 

X  JL  X 

^  o,  and  we  take  it  for  granted  that  the  integral  equation  then  has  a  unique 
solution  which  is  positive  everywhere o 


15 


For  simplicity  write  the  integral  equation  (Io7)  in  the  form 


Uo8) 


(p  =  Acp  +  B  , 


where  A  is  a  linear  integral  operator  and  B  is  a  known  function 0  Let  cp  be 
the  solution  of  (Ic8)  and  let  cp^  =  (p  +  ^  be  a  zeroth  approximation 0  Neither 
the  flux  (p  nor  the  error  is  known*  only  their  sum  (pQ  is  given*  Then  the 
successive  approximations  are  defined  by 


(1.9) 


fn  =  %-l  +  B 


The  error  is  itself  the  solution  of  an  integral  equation  like  (lo8),  for 
we  have 


(1-10) 


€0  =  A€0  +  ^0  -  ?! 


Bounds  for  the  error  e Q  can  be  estimated  from  (l*10)o 

The  errors  in  the  successive  approximations  are  obtained  as  follows: 


(loll) 


=  cp1  -  cp  =  A(cp0  -  <p) 


A€1  =  A  €o 


Hence  the  error  — 0  if  and  only  if  A11^  -~0o  The  important  fact 
is  that  if  it  is  possible  to  derive  an  estimate  |«J  -  ^'€nJ  '  0<r)  Kl> 
from  the  equality  then  it  can  be  asserted  that  the  prror  is  reduced 

by  at  least  the  fraction  I-77  at  each  iteration  Together  with  an  estimate  of 
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INTERNAL  REGION,  R2 


FIGURE  2 

SINGLE  CONVEX  INTERNAL  REGION 


€  ,  one  is  then  in  a  position  to  determine  a  bound  for  the  error  in  any  iterate. 
Clearly  the  convergence  of  the  successive  approximations  to  (p  is  equivalent 

to  £  —  0. 

n 

4 o  Convex  Regions 

This  section,  and  section  7  below*  are  concerned  with  the  problem  of 
convergence  of  the  successive  approximations  defined  by  (l09)0  For  a  single 
convex  internal  region  tne  question  of  convergence  is  answered  by  the  theorem 
on  page  22  0  The  present  section  contains  the  proof  of  this  theorem*  as  well  as 
derivations  of  inequalities  giving  bounds  for  the  flux  and  the  errors o  Non- 
convex  regions,  for  which  less  is  known,  are  discussed  in  section  7 <,  By  way  of 
examples,  the  flux  distribution  in  a  slab  water  hole  is  calculated  in  section  5, 
and  in  a.  cylindrical  water  hole  in  section  6* 

A  convex  region  is  defined  by  the  property  that  every  chord  intersects 
the  surface  at  most  twice o 

Let  Rg,  Figure  2,  be  convex.  The  flux  in  is  the  solution  of 

(1,12)  iy>  =  (DrD2)  f  (v«i)G 

^2  Pl2  x+2  R2 

With  the  physical  situation  in  mind  we  introduce  the  following  as¬ 
sumptions  and  notations^ 


(DrD2)/D2  =  a  >  0  , 


2  2 

1  2 


p  ^  0 


Other 

cases  can  be  discussed  in  a  manner  similar  to  the  analysis  below.  Also  observe 
that  G>0,  dG(x)/dx  <0,  and 


r  /, 


Sri  + 


(qo-qJQ  /  >o 
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The  integral  equation  (lo!2)  now  reads 


(I- 13) 


9  =  a  9r+P  /  9^  +  Y 


.  To  apply  the  method  of  successive  approximations  we  require  a  starting 
function  cp^  and  an  estimate  of  the  error,  -  cpQ  -  cp0  These  are  most  easily 
obtained  by  determining  upper  and  lower  bounds  for  the  flux  (p0  The  simplest 
bounds  are  the  largest  and  smallest  asymptotic  values  of  the  flux  but  better 
bounds  can  usually  be  calcuiatedo  Consider  region  R0  and  let  (p  „  and  (p  < 
denote  the  maximum  and  the  minimum  values  of  <p  there o  The  convexity  assumption 
implies  0G/0n  <0  on  S^o  Since  a,  p,  and  y  are  non-negative,  it  can  be  con¬ 
cluded  from  (lol5)  that 


(iol4) 


<P  <  CKp, 


fm  + 

Jrl-n 


min  J  5n 
S2 


G  +  y 


(Icl5) 


(p  >  ocp, 


max  J  0n 


+  P9r 


G  +  y 


Inequalities  (Id4)  and  (lol5)  furnish  upper  and  lower  bounds  for  (p 

at  any  point  r  in  region  R0  in  terms  of  q)  and  cp  .  there  o  In  particular 
“  fc  Tmax  Tmin 

let  r  be  a  point  at  which  <Pmax  is  attainedo  Solving  (io 14)  and  (l015)  for 


cp  gives 
Ymax 


G  +  y 


(I«16) 


./ 


^min  J„  dn  +  ^ 

SZ 


1  -  p 
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Similarly , 


(1.17) 


bounds  for  cp  .  are 
Tmin 

[  5G 

a<fynex  I  9  n  ^ 

S2  _ 

1  "  P  /  G 

R0 


W>, 


max 


-  <p  .  - 

Ymin 


( 


G  +  y 


1  -  a 


f  QG 

L  5n 


o 


The  inequalities  (lol6)  and  (1*17)  may  be  manipulated  in  various  ways, 

The  extreme  inequalities  can  be  solved  to  provide  a  lower  bound  for  <Pmax  snc^  an 

upper  bound  for  <p  .  o  Since  it  is  apparent  that  <p  .  -  min(q/£),  an  upper 

min  *  min 

bound,  w,  for  (p  can  be  obtained  from  (lol6)„  Then  an  improved  lover  bound, Y, 
T  7  Tmax 

for  ©  .  is  riven  by  the  left  side  of  (l<,17)  „  To  what  extent  further  manipula- 
^min 

tion  is  desirable  depends  on  the  problem  at  hand* 

It  should  be  observed  that  the  upper  bound  for  (p  -n  furnished  by  (l„16) 
is  itself  a  significant  number.  If  q<,/S  ^  >  q^/2^  then  the 

facts  that  Dg  -  and  that  Rg  is  convex  imply  <pmin  occurs  on  Also  since 

R  is  convex  and  is  large  enough  for  diffusion  theory  to  apply,  the  flux  would 
not  be  expected  to  vary  greatly  on  the  surface  S^o  These  ”facts,!  imply  that  the 
upper  bound  for  (pmin  already  furnishes  an  estimate  of  how  much  flux  peaking  will 
occur  in  the  external  region » 

As  a  zeroth  approximation  to  the  flux  we  can  take  (p^  =  (ty  +^)/2  and 
know  the  error  ^  does  not  exceed  -  ^)/2o  The  error  in  (p^  can  be  esti¬ 
mated  in  various  ways.,  depending  on  how  many  iterates  are  available ,  Operating 
on  (lolO)  with  A  gives  ^  =  k*1  +  <PX  -  from  which  it  is  possible  to  estimate 
€  much  as  (p  is  estimated  by  (1*14)  and  (lu15)o  Alternatively,  from  cpx  +  ^ 

=  Vq  +  €q  are  derived  the  inequalities 

min  €0  +  <P0  -  ^  -  ei  "  “a*  e0  +  90  “  ° 
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Perhaps  the  simplest  and  most  useful  bounds  for  are  obtained  from 


a 


/ 


€  9G  +  o 
0  3n  ^ 


9 


which  furnishes 


<  -  a  max 


+  p  max 


9 


or,  upon  rearrangement, 


o 


The  error  will  tend  to  zero  and  the  successive  approximations  will  converge 
provided 


(Ida) 


(P  -  a*-2)  f 

h 


G  +  a  <  1 


o 


It  is  instructive  to  solve  (l„18)  for  ae 


The  preceding  inequality  is  certainly  satisfied  if  a  <  1* 
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So  far  we  have  proved  the  first  part  of  the  following  theorems 


Let  a,  jjj,  and  %  be  non-negative  and  let  the  region  R^  be  convex,  then 
if  a  <  1,  the  integral  equation  (idS)  can  be  solved  by  the  method  of 
successive  a p proximat ions 0  But  if  a  >  1,  the  successive  approximations 
can  fail  to  converge  for  all  sufficiently  small  regions  R^j  or  for  a 


given 


Rg  the  approximations  can  fail  to  converge  for  all  suffix 


ficiently  large  a. 


Before  proving  the  latter  assertions  we  shall  illustrate  them  in  a 
particular  case0  Suppose  p  =  0  so  that  we  have  merely 


. .  /. 

n  -1  n- 


3G 
■1  3n 


Now  let  €  -  S  >  Oo  Then  s  <  a8 

n-l  n 


=  a  S  (-1  +  j  G  )  o  Choose 


a  >  1  and  Rp  so  that  €  <  -  B0  This  will  certainly  be  true  if 

n 

J  G  <  (a  -  l)/a„  But  then  ^n+1  t  -a  8  J  ^  S  y  and  the  successive 
R2  S2 
errors  will  never  tend  to  zeroj  they  will  in  fact  increase „  It  seems  remarkable 
that  an  arbitrarily  good  approximation,  i0e«,  8  as  small  as  you  please,  cannot 

ensure  the  convergence  of  the  iterates □ 

A  slight  elaboration  of  the  argument  just  given  will  complete  the 
proof  of  the  theorem 0  Let  ^  satisfy 

0  <  81  s  Vl  5  S2  * 


M  =  max  J  G  , 

h 


m  =  min  /  G 


Then  vre  have 


-a  +  (aKj2  +  pS1)m  ^  eR  <  -aS1  +  (a/^8.^  +  p  82)M 


and  if  M  and  a  satisfy 


M  < 


(a  -  1)  8 


*W8*  ’ 


we  find 


-8*  S‘.  S-»l 


Now  we  want  to  show  that  e  ^  and  that  the  argument  can  be 

continued o  Clearly 


*Z  \  -  (  \K^  +  P  S8)m7  -  €1  ■’  -  °2 


Since 


<  Vi‘  -  e  VM  s  ^S1 


we  have  finally 


8!  5  Vi  £  -*  8e 


Hence 


^2  <  *n+2  “  ^1 

S1  "  €n+3  "  “4  8  2 


The  errors  remain  outside  the  interval  (-  8),  and  the  iterates  do  not  converge 

This  completes  the  proof  of  the  theorem,. 
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2k 


SLAB  WATER  HOLE  IN  AN  INFINITE  REACTOR 


Now  that  the  discussion  of  successive  approximations  for  a  convex 
internal  region  is  complete,  we  illustrate  the  application  of  the  integral 
equation  by  determining  the  flux  distributions  in  slab  and  cylindrical  water 
holes o 

5 o  The  Slab  Water  Hole  in  an  Infinite  Reactor 

In  all  our  examples  we  shall  use  the  following  physical  constants: 


Reacting  Material,  R~^ 

Water  Hole,  R^ 

D 

0275 

.162 

2 

K 

.200 

.121 

q 

1.000 

2.009 

Since  a  =  (D^  -  =  .70,  the  successive  approximations  are  sure  to  convergee 

Let  the  water  hole,  R^,  be  a  slab  4  cm  in  thickness  and  let  the  reac¬ 
ting  material,  R^,  extend  to  infinity  as  shown  in  the  sketch  on  page  24 0  The 
flux,  determined  by  an  analytic  solution  of  the  differential  equation,  is 

f  44.7  e-“447x  +  18.2  ,  x  >  2 

(I  >19)  cp  =  j 

-53ol  cosh(.347x)  +  102.9  ,  |xl  ^  2  . 


This  is  illustrated  in  Figure  3,  page  26« 

The  Green 3 s  function  for  the  problem  is 

lx-x3l 

G(x,x«)  =  e  /2\ 


and  the  integral  equation,  for'  (p  in  R^, 


reads  after  slight  modification: 


i 


t 

i 
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(1.20)  <p(x) 


cosh^x  + 


«p(2)e 


-2*1 


+ 


2K1D2 


x-x 


dx»  + 


D  * 
^2  1 


One  notes  that  the  surface  integral  can  be  explicitly  evaluated  in  those  cases 
in  which  cp  is  constant  on  the  surface. 

We  shall  calculate  bounds  for  the  flux  and  then  obtain  an  approximate 
solution  of  (lo20)  by  iteration.  As  was  pointed  out  on  page  20,  the  maximum  and 
minimum  values  of  the  flux  in  the  water  hole  occur  at  x  =  0  and  x  =  2.  The 
right  side  of  (1.20)  is  made  smaller  if  (p(x)  is  replaced  by  cp(2) .  Then  putting 
x  =  2,  (1.20)  can  be  rearranged  to  provide  the  lower  bound  cp(2)  -  35.4,  Next 
put  x  =  0  and  replace  (p(2)  by  the  lower  bound  35. 4o  This  gives  <p(0)  —  51. 3o 
An  upper  bound,  37.2,  for  (p(2)  is  obtained  by  putting  x  =  2  and  replacing 
cp(x')  by  its  majorant,  51.3.  Finally  cp(0)  £  47.6  follows  on  replacing  <p(2>  by 
37.2  and  (p(x9)  by  35.4.  In  summary 

47.6  <  cp(0)  <  51.3  f  Cp(0)  =  49.3  | 

35.4  <  <p(2)  <  37.2  1  <p(2)  =  36.5  J 


It  is  in  a  way  surprising  that  such  narrow  bounds  for  the  flux  at  the 
center  and  at  the  boundary  are  obtained.  This  fact  results,  perhaps,  from  the 
geometry  which  tells  us  that  the  flux  has  a  constant  minimum  value  at  the  sur¬ 
face  of  the  water  hole,  -tor  many  purposes,  the  fact  that  (p  on  the  surface  lies 
between  35  and  37  is  all  the  information  required.  Since  the  asymptotic  value 
of  the  flux  is  18.2,  the  peaking  effect  of  the  water  hole  is  well-defined. 

A  first  approximation,  cp^,  to  the  flux  has  been  calculated  in  two  ways 

First,  put  cpQ(x)  =  43.4,  -2  $  x  ^2,  and  (pQ(2)  =  36.  The  values  of  (p^x)  cal¬ 
culated  from  (1.20)  are  shown  as  circles  in  Figure  3  (page  26).  The  approxima¬ 
tion  (p  =  43.4  is  very  rough,  yet  (p^(x)  turns  out  to  be  close  to  the  actual  flux 
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REACTOR  MATERIAL 


CYLINDRICAL  WATER  HOLE  IN  AN  INFINITE  REACTOR 


This  is  due  in  pciri  bo  bhe  fact  "that,  a  good  approximation  lor  the  llux  at  the 
surface  is  available 0 

A  much  better  initial  approximation,  <p  ,  should  be  given  by  a  line 
segment  passing  through  the  flux  values  50  and  36  at  x  =  0  and  x  =  2  respec- 
tively,  For  cpQ(x)  we  have  <pQ  =  -  7 1  x  1  +50,  with  Ixl  <  2,  The  values  of 
<p^(x)  obtained  in  this  way  are  shown  as  squares  in  Figure  3  (page  26) „  Evi¬ 
dently  the  more  elaborate  estimate  of  (p^  was  not  worth  the  ef forte 

6  o  The  Cylindrical  Water  Hole  in  an  Infinite  Reacaor 

The  flux  distributi.on  in  a  cylindrical  water  hole  (see  sketch  on  page 
28),  8  cm  in  diameter,  was  examined  by  means  of  the  integral  equation,,  The  Green's 
function  is  r  -  r'|  )/2tt;  the  constants  are  the  same  as  in  the  slab  treated 

above o  For  the  integral  equation  we  now  have 


(I, 21) 


<p(r) 


cpCr'jKQ^r 


r«|) 


27rD, 


R, 


Kq^Ii;  -  r'|  )  + 


Vl 


o 


The  surface  integral  in  (I, 21)  is  transformed  to  a  volume  integral 
by  Green's  theorem,  Expansion  of  K^(kJ  r  -  r9|  )  by  the  addition  theorem  for 
Bessel's  functions  and  integration  over  the  angular  variable  replaces  the  double 
integrals  by  simple  integrals  which  are  manageable. 

The  exact  analytical  solution  of  (I <>21)  is 

r  -45o5I  (o347r)  +  102,9  ,  r  <  4 

(I»22)  <p(r)  =  \ 

98,3  K0(„447r)  +  18,2  ,  r  >  4 
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Neutron  Flux  Level 


Bounds  for  (pmin  =  (p(4)  and  (pmax  =  <p(0)  are  found  exactly  as  before. 


51,6  <  <p(4)  ^  54o5  , 
52  o  2  ^  (p(0)  ^  60.5  , 


(p(4)  =  32.8 ' 

ffi(0)  =  57.5  J 


The  first  approximation  (p^(r)  was  calculated  from  cpQ  —  46,  cp(4)  —  33.  The 
iterated  flux  cp  ,  as  well  as  the  correct  flux  <p,  is  illustrated  in  Figure  4 
(page  30) . 


7.  Non-Convex  Regions 

The  question  of  convergence  of  the  successive  approximations  is  by 
no  means  simple  for  regions  that  are  not  convexo  For  purposes  of  computation 
it  is  not  only  the  fact  of  convergence  which  is  important}  equally  important  is 
the  requirement  that  the  successive  approximations  converge  rapidly.  It  is 
easy  to  see  that  no  general  criteria  can  ensure  the  convergence  in  a  manner 
which  is  useful.  In  equation  (1.13),  let  fi  =  0.  We  prove: 

For  every  a  >  0  there  exist  regions,  Rg,  such  that  the  iterated  error 
is  at  some  point  arbitrarily  large .  although  e_Q  is  as  small  as  you  please. 


There  is  no  loss  in  generality  in  allowing  to  be  discontinuous. 
With  respect  to  some  point  in  a  region,  Rg,  let  €q  =  8  >  0  on  that  part  of 
the  surface  where  3G/3n  >  0  and  let  =  -8  where  3G/3n  <  0.  Then 


Clearly  there  exist  non-convex  regions  such  that 


is,  for  some  point. 


greater  than  any  preassigned  number.  This  proves  the  assertion.  Of  course, 
the  counter  example  is  pathological}  one  can  hope  that  ordinarily  the  process  of 
iteration  will  not  lead  to  disastrous  results. 
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FIGURE  5 

CROSS-SHAPED  WATER  HOLE  IN  AN  INFINITE  REACTOR 


What  we  can  say  is 


that  €  will  certainly  tend  steadily  to  aero  if 
n 


(1-25) 


max  € 


*  «  /  If  I  ♦  >  / 


G  <  1 


The  only  way  of  improving  on  (1.23)  is  to  obtain  some  knowledge  of  the  distri¬ 
bution  of  the  error  on  the  surface .  One  way  of  doing  this  is  to  obtain  bounds 
for  cp  at  various  points  on  the  surface  from  the  inequalities  analogous  to  (1.14) 
and  (Iol5)  which  are  applicable  to  non-convex  bodies.  These  inequalities  will 
be  derived,  for  the  cross-shaped  water  hole,  in  the  next  section. 

8,  The  Cross-Shaped  Water  Hole  in  an  Infinite  Reactor 

Our  interest  here  is  to  show  one  way  in  which  bounds  for  the  flux  in 
non-convex  regions  can  be  calculated.  The  cross-shaped  water  hole  is  used  only 
for  purposes  of  illustration.  We  already  know  an  upper  and  a  lower  bourn  for 
the  flux— the  asymptotic  flux  values.  Other  bounds  can  be  obtained  by  modifying 
the  physical  constants  to  yield  simpler  integral  equations  whose  solutions  are 
everywhere  greater  than  or  less  than  the  correct  flux.  Our  procedure  here  will 
be  to  obtain  inequalities  analogous  to  those  derived  in  section  4  for  convex 
regions.  Let  the  point  of  interest  be  r.  Figure  5  (page  32),  and  draw  the 
auxiliary  surfaces  shown  by  the  dotted  lines.  In  addition  to  the  cross  and 

surface  Sg,  we  consider  the  six  regions  Rj,  ...  ^  and  surfaces  Sj,  ...  S^. 
All  normals  are  directed  outward o  From  the  equality. 


!2+SIII+SV 


sii+siv 


SII+SIII+SIV+SV 


follow 
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(lo24)  <PT-  <  cp  . 


5n  +  ^min 


m_(D 

6n  vmax 


S2+SIII+SV 


SII+SIV 


sii+siii+siv+sv 


(1.25) 


-  % 


On  unax 


On  uni n 


S2+SIII+SV 


SII+SIV 


SII+SIII+SIV+SV 


where,  in  the  first  instance,  the  maximum  and  minimum  values  of  (p  are  to  he 
taken  on  the  appropriate  surfaces 0 

These  inequalities  appear  to  better  advantage  when  rearranged  thus 


(lo26) 


OG  .  /  0G  .  ,  . 

^  On  ~  ^min  I  On  ^min  ^max^ 


SII+SIII+SIV+; 


(1.27) 


(p  —  >  (p  /  *—  +  ((p  -  (p  ,  ) 

Y  On  unax  /  On  v  unax  unxn7 


SII+SIII+SIV+SV 


From  (lo26),  (Io£7),  and  the  integral  equation  (lo!3)  are  obtained  the 
bounds,  for  <p  in  the  cross. 


(Xo28)  (p  -  acp  .  /  ~  +  a(cp  .  -  cp  ) 

t  unm  /  On  'unxn  unax7 


^  +  p<p  /  G  +  y 
On  Kunax  /  r 


SII+SIII+SIV+SV 


(Io29)  cp  >  acp 


max  /  On 


+  a(<p  -  cp  .  ) 
unax  unm' 


SII+SIII+SIV+SV' 


to  +  p<pmin  J  G  +  Y 
*2 
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CROSS-SHAPED  WATER  HOLE  IN  AN  INFINITE  REACTOR 


The  right  hand  sides  of  (lo28)  and  (I.29)  are  unambiguous  if  <Pmax  and  ^min  are 
understood  to  be  the  maximum  and  minimum  values  of  (p  in  the  entire  region 
Rj  +  «o*o  +  the  inequalities  are  sharper  if  (p^.^  and  (pmay  are  interpreted 

to  best  advantage  each  place  they  occur. 


No w  we  shall  suppose  that  cp  occurs  at  the  center  of  the  cross,  r, , 

Tmax  ”1 

and  that  cp^.^  (in  the  cross)  occurs  at  the  corner  r^,  Figure  6  (page  35).  At 
the  point  r^,  6G/0n  is  everywhere  negative,  i.eo,  the  cross  looks  convex  from 
the  center.  Therefore  inequalities  (I016)  apply  with  cpmax  =  <p(r^)  and  = 

The  only  lower  bound  available  for  <p  in  regions  R_^  and  is  the 
asymptotic  value  q^/2  The  sharpest  bounds  for  (p(iv,)  are  obtained  by  using 
(I„24)  and  (1.25): 


(lo31)  q>(r2)  ^  ®pC£q) 


+  p<p( 


^  / 


G  +  y 


S2+SIII+SIV+SV 


SIII+SIV' 


Solving  (1.30)  and  (l„31)  for  <p(r2)  then  gives  the  inequalities  we  set  out 
to  derive; 


With  these  cumbersome  inequalities  we  end  our  discussion  of  non-convex 
regions  o  In  a  later  report  we  hope  to  give  comprehensive  numerical  results  il¬ 
lustrating  the  methods  discussed  in  this  report  in  their  application  to  cross¬ 
shaped  water  holes □ 

Before  turning  to  the  very  important  method  of  superposition,  we  make 
an  observation  on  the  boundary  values  satisfied  by  the  flux  in  the  internal 
region  R^o 

90  Appendix  to  Part  I 

As  was  pointed  out  before,  the  integral  equation  (lol£)  for  the  flux 
in  R^  can  be  solved  without  any  consideration  of  the  flux  in  R^  Since  <p 
satisfies  the  differential  equation  (Io3)  it  is  natural  to  ask  whether  or  not 
the  integral  equation  (Id2)  is  equivalent  to  the  differential  equation  and  a 
boundary  condition  on  which  involves  only  the  flux  in  R^o  This  question 
will  be  examined  here® 
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From  the  equations 


(1.3) 

DgV2<p  -  Xg(p  +  q. 

(1-5) 

D-j^V^G  -  X-jG  : 

is  obtained  at  once  the  integral  equation 

(1.33)  D2,p  =  d2  /<G  f  -  <p  f )  +  (*/-  fat  +  q2  fa  . 

ft  *2  h 


Certainly  this  integral  equation  cannot  determine  (p  since  no  boundary  condition 
has  been  imposed,.  Now  equation  (I „ 33)  could  be  transformed  to  yield  (l„12)  if 
the  boundary  condition  were  known..  Hence  the  boundary  condition  must  be  im¬ 
plicit  in  the  equation  obtained  by  subtracting  (l„33)  from  (l„12).  There 
results 


=  0 


or 


(1=34) 


The  preceding  line  is  indeed  a  boundary  condition,  a  very  complicated  boundary 
condition,.  It  asserts  that  we  must  select  from  all  the  solutions  of  (1  =  3)  the 
one  whose  boundary  values  make  (l„34)  an  identity  as  a  function  of  r. 
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Suppose  S„  and  S  are  concentric  cylinders  with  radii  a  and  b  respec- 
tivelyo  Then  <p  and  9(p/3n  are  constant  on  and  the  boundary  condition  reads 


/V<a>  -  tt-7  M  -  /°  »  1  -  /G 

1  S2  S2  "1+2 


Now  put  r  =  Oo  Then 


i  -  Kn(Klb)  v 

G(r')  =  2~rr  L  K0(Klr,)  "  Ig^b)  W')-7  * 


„  _  K(Kb) 

I  -  -  *A<V>  -  i^5T  Ii(Kia)-7 


so  that  the  boundary  condition  becomes 


(l<>35)  Dxk 


K_  (k,  a)  I-,(K-,a)  q. 


K0(«i«)  W> 


+  V  K^bT  -  I^b)  -7«"(a)  -  -  I^bT 


In  particular,  if  b  =  oo  , 


^1 

Dl*lKl(*la)  +  Ki^Kia)  =  0 


Of  course,  the  boundary  condition  (Io35)  may  be  derived  in  a  completely 
elementary  manner  0 
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II o  THE  METHOD  OF  SUPERPOSITION 


Numerical  solution  of  the  integral  equation  for  complicated  configura¬ 
tions  is  a  formidable  tasko  The  difficulties  are  of  two  kinds 0  The  first  of 
these  is  typified  by  the  problem  of  many  water  holes  in  an  infinite  reactor,  or 
even  of  a  single  water  hole  with  a  complicated  shape 0  Difficulties  of  another 
kind  enter  into  the  calculations  when  the  reactor  does  not  extend  to  infinity* 
For  then  we  are  faced  with  an  unwieldy  Green 3 s  function,  if  it  is  known  at  all, 
and  with  source  terms  which  cannot  be  considered  constant  in  each  homogeneous 
part  of  the  system*  These  obstacles  are  largely  overcome  by  the  use  of  the 
method  of  superposition,  a  procedure  which  has  often  been  applied,  implicitly 
if  not  explicitly* 

Consider  two  water  holes  in  an  infinite  reactor*  Let  (p^  and  <p^ 
be  the  flux  distributions  due  to  each  in  the  absence  of  the  other  one*  Let 
qf/^i  be  the  asymptotic  value  of  the  flux*  By  superposition  we  mean  the  flux 
distribution 


<p(s)  =  ♦  <p<3>  -  <qi/ 21)  . 

(s) 

The  question  to  be  answered  is  how  well  does  <pv  7  approximate  the  correct  flux 
distribution  <p* 

The  term  superposition  will  be  used  also  in  another  sense*  Let  be  ■ 
the  flux  distribution  for  a  water  hole  in  a  reactor  which  extends  to  infinity* 
Can  the  flux  distribution  for  the  same  water  hole  in  a  finite  reactor  be 
reasonably  well  described  by  the  product 

<p^  - 

where  F  is  a  suitable  weighting  function,  usually  easily  determined?  We  shall 
refer  to  the  two  kinds  of  superposition  as  superposition  by  addition  and  super¬ 
position  by  multiplication* 
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REACTOR  MATERIAL 


APPROXIMATION  TO  THE  CROSS-SHAPED  WATER 
HOLE  BY  CYLINDERS 


More  convincing  than  any  abstract  mathematical  discussion  of  the  power 
of  superposition  is  the  drastic  test  provided  by  breaking  up  the  cross  shaped 
water  hole  into  cylinders  and  superimposing  the  flux  distribution  of  the 
cylinders o  The  cross  section  of  the  cross  shaped  water  hole  we  are  considering 
can  be  thought  of  as  made  up  of  13  squares.  Figure  7«  The  method  of  superposi¬ 
tion  used  was  to  replace  each  square  by  a  circle  of  equal  area  and  to  superim¬ 
pose  the  flux  distributions  of  the  circles .  The  correct  flux  values  as  well  as 
the  values  obtained  by  superposition  are  plotted  in  Figure  8  (page  43) 0  In  this 
case  superposition  underestimates  the  ilux0  The  error  in  (p  -  (q^/2^),  the  flux 
rise  over  the  asymptotic  value,  is  less  than  20% 0  As  will  be  shown  below,  the 
fact  of  underestimation  can  be  predicted  and  a  correction  can  be  appliedo 

Formulae  which  justify  the  method  of  superposition  and  which  can  be 
used  to  estimate  the  superposition  error  will  now  be  derived 0  Also  a  simple 
first  order  correction  to  superposition  is  obtainedo  In  section  1,  superposi¬ 
tion  by  addition  is  discussed  and  illustrated  by  finite  slab  water  holes  in  an 
infinite  reactorQ  Multiplicative  superposition  is  discussed  in  section  2,  where 
we  take  as  our  example  a  slab  water  hole  in  a  finite  slab  reactor o  In  section  3 
a  finite  reactor  with  internal  and  external  reflectors  is  considered  briefly 0 

la  Superposition  by  Addition 

We  shall  examine  the  flux  distributions  created  by  two  diffusing 
regions  and  in  an  infinite  medium  R^,  as  illustrated 0 
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NEUTRON  FLUX  DISTRIBUTION  AT  A 
CROSS  SHAPED  WATER  HOLE 
IN  AN  INFINITE  REACTOR 


(Ilol) 


where 


The  flux  <p  is  the  solution  of  the  integral  equation 


iyp  =  u<p)  +  q 


di.2)  l(»)  =  (h.rz)  /q>  f  +  -  d3)  Af 

S2  S3 

+  (''i2  -  Jv*  +  /a P  > 


«  =  /q-,0  +  j  (%  -  q,)G  +  /(q,  -  q,)a 

W  R3 

(g)  (3) 

Let  cp  and  <p'  '  be  the  flux  distributions  due  to  media  and 
Then 

=  (Di  -  op  f<i p(j)  +  (k*  -  K*)n.  /<*<*> 


/qxG  +  /(q,  -  q,)G  , 

w3  Rj  j 


i  =  2,3 
i  =  1,  or  j 


Let  cp^  -  +  <p^  -  q1/^1  be  the  superposed  flux.  We  want  to  estimate 

8<p  =  <p  -  <p's\  It  can  be  verified  that  S  cp  satisfies  the  following  integral 


equation 


(IIo3) 


D.  8<p  =  L(Scp)  +  8q  , 


i  =  1,2,3  , 
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where  the  operator  L  is  defined  in  (H°2)  and 


Equation  (II. 3)  is  an  inhomogeneous  integral  equation  which  differs 
from  (II.1)  only  in  the  inhomogeneous  term  8q.  If  8  Q  were  zero  then  8<p  =  0 
and  superposition  would  give  the  correct  answer.  Note  that  (II.4)  states  that 
the  extent  to  which  8Q  differs  from  zero  depends  on  the  deviations  of  the  fluxes 
<P(2)  and  <p^  from  their  asymptotic  values  in  regions  Rg  and  Rg  respectively. 

The  superposition  flux  <p^  is  that  approximation  to  (p  in  which  8<p  =  0. 
Now  8<p  =  0  can  be  thought  of  as  a  zeroth  approximation  in  solving  (II. 3)  by 
iteration;  hence  a  first  approximation,  Scpj,  to  8<p  is  obtained  by  putting 
8  <p  =  0  in  the  right  hand  side  of  (II. 3).  There  results 

89X  =  Bq/d.  > 

but  it  does  not  necessarily  follow  that  <p^  +  8  is  a  better  approximation 
to  <p  than  <p^ .  It  can  only  be  hoped  that  the  iterate  is  more  nearly  correct 
than  thb  zeroth  approximation o  The  correction  is  in  the  right  direction,  for 
it  can  be  assumed  that  Sep  will  have  the  sign  of  8q,  at  least  if  8q  is  either 
everywhere  positive  or  everywhere  negative. 

The  sign  of  8q  in  the  internal  regions  can  be  determined  for  certain 
problems  of  importance.  Let  the  point  of  interest,  r,  be  in  Rgo  Now  we  make 


45 


the  assumption  Rg  is  convex  so  that  5G/9n  is  negative  on  Sg .  Put  (p ^ ^ -  (q^/Z 
=  Acp^,  Tnen,  if  a  =  (D-^  -  Dg)/Dg  -  0,  p  =  K-^  -  -  0,  we  have 

^  <  a  /  max^A  (p(3)  -  ming  A  J 

+  /~P  max^  A(p(5)  +  K^a  ming  A  (p^5^_7  J  G  +  —■  , 


>  a  /"mii^  Aq/3)  -  maXg  A 

+  jT$  min^  A(p^3^  +  K-^ct  maXg  A  (p^_7  ^G  +  ■—  , 


where  €  is  the  contribution  of  the  terms  containing  the  integrals  over  and 
Rgo  In  deriving  these  inequalities  it  must  be  remembered  that  our  assumptions 
imply  A  (p(3)  >  0. 

For  a  >  0  it  follows  that  8  Q  <  0  in  if 


a  (3)  1  -  k  Z  J  G 

maxR  A  <p^  1  J 

(II»5)  - -a-t?y  - o’  1  ~f  ~ 

ming  A  <p'  '  1  +  *  G 


ming  A  (p^S^  (a  +  p  J  G) 
2  FL 


and  S  Q  -  0  in  Rg  if 


(II. 6) 


rnaxg  A  fW  1  +  |  [  G 


A(p(5ha  +  p  f  G) 

3  </d 
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The  criteria  (II, 5)  and  (II, 6)  are  exact,  but  their  usefulness  de¬ 
pends  upon  an  estimation  of  <=.  In  many  cases  the  integrals  over  S3  and  Rg  will 
be  negligible  compared  with  those  over  Sg  and  Rg  so  that  the  terms  containing  € 
can  be  ignored.  The  important  contributions  to  the  integrals  over  S,  and  R„ 
will  come  from  that  region  of  R^  nearest  to  R^  since  G,  9G/0n,  and  Acp  all 
decrease  rapidly  with  distance 0  Thus  an  approximation  to  €  is  obtained  if  one 
replaces,  in  both  integrals,  A<p^  by  a  suitable  average  which  is  ap- 

proximately  the  value  at  the  nearest  point o  Then  we  find 


( I,  -  fa 


and  e  has  the  sign  of  2-^  -  2^. 

The  inequality  (II. 5)  can  be  satisfied  only  if  e  is  negative  and,  in 
absolute  value,  not  too  small.  On  the  other  hand,  it  may  be  possible  to  satisfy 

(11.6)  when  the  €  term  is  discarded.  If  €  >0,  it  can  be  discarded,  and  (ll,6) 
then  provides  a  sufficient  condition,  easily  evaluated,  in  order  that  S  Q  -  0, 

On  physical  grounds,  and  according  to  the  examples  which  have  been  calculated, 
the  assumptions  on  which  (11,5)  and  (lIo6)  are  based,  namely  a  >  0  and  (5  -  0, 
imply  that  8(p  and  hence  8  Q  are  greater  than  zero.  This  means  that  inequality 

(11.6)  rather  than  (II. 5)  should  apply.  Also,  since  >  0,  we  have 

€  >  0o  In  other  words  (Ho5)  will  probably  never  be  satisfied  and  the  suffi¬ 
cient  condition  (11,6)  may  be  satisfied.  Similarly,  if  a  and  p  are  negative  so 
that  the  flux  is  depressed  in  Rg  only  one  of  the  ensuing  inequalities  will  in 
fact  be  applicable. 

It  is  of  some  interest  to  compare  the  superposition  approximation  with 
the  analytic  solution  in  certain  simple  cases  where  the  analytic  solution  may 
be  obtained. 
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WATER  I  Y WATER 


TWO  SLAB  WATER  HOLES  IN  AN  INFINITE  REACTOR 


(a)  Two  slab  water  holes  in  an  infinite  reactor 

We  first  examine  superposition  for  the  case  of  two  slab  water  holes > 

4  cm  thick,  and  4  cm  apart  in  an  infinite  reactor  (see  sketch  on  page  48).  The 
correct  flux,  (p,  and  the  superposition  flux  illustrated  in  Figure  9  (see 

page  50)  are  given  respectively  by  the  formulae! 


(p  =  18.2  +  15„8  cosh  .447x  , 

=  102.5  -  6.62  eo548x  -  97.1  e~o548x 


0  <  x  <  2 


2  <  x  <  6 


=  18.2  +  276.9  e 


-o447x 


x  >  6 


=  18.2  +  14.9  cosh  .447x 


0  <  x  <  2 


=  102.5  -  52.8  cosh/7348(x-4^7  +  44.6  e“°477(x+4)  2  <  x  <  6 


=  18.2  +  274.0  e 


-»447x 


x  >  6 


(b)  The  reflected  slab  reactor 

We  next  examine  superposition  for  a  finite  reflected  slab  reactor. 
Figure  10.  The  flux  <p,  the  solution  of  Dtp"  -  2cp  +  q  =  0,  is  given  by  the 


formulae 


(p  =  c^  cosh  k^x  +  ’ 

=  c2  e-^lxl+  q2/£2  , 


D2/<2^q2/,^2  "  ql/^l^ 


C1  Dg^cosh  k^&  +  D^x^sinh  ^a 


c2  =  -  e*2a  sinh  Kla 


x  in  L 


±  This  example  was  suggested  by  Mo  Danzker,  Reactor  Theory  Section,  Physics 
Department,  Westinghouse  Atomic  Power  Division 


Neutron  Flux  Level 


FIGURE  9 
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The  superimposed,  flux,  obtained  by  combining  the  fluxes  when  first 
one,  then  the  other  reflector  is  replaced  by  the  reactor  material,  is  given 


jp(s)  _  la  cosh  k^x  +  ^  > 

=  a,e-«l(^)+  + , JX  , 


in  R^ 


x  in  IV 


with 


di  = 


d2  = 


“  ql/^l) 


\\  +  hKl 


Dl*l  ^2^  2  2  ~  qj/^  1^ 


°2*2  +  h*l 


The  ratio  (<p  -  cp^S^)/((p  -  q^/2^)  has  been  plotted  in  Figure  10, 
page  52 ,  for  several  sets  of  values  of  the  physical  constants »  In  region 
this  ratio  is  the  constant , 

i  -  (j2vvp  ^ 
i  +  a>2yi>1«1) 

so  that  superposition  gives  the  correct  answer,  underestimates  or  overestimates 
the  flux  in  the  core  according  as  =  ^  ^  ^  or  ^  simple  rela~ 

tionship  holds  in  the  reflectorG  Indeed,  as  shown  by  curve  V,  Figure  10, 
cp  -  <p^  can  change  sign0 


£o  Superposition  by  Multiplication 

Consider  a  reactor  of  finite  extent  composed  of 
two  homogeneous  parts,  R^  and  R^o  We  have  in  mind  a  bare 
reactor  with  a  water  hole0  The  reflected  reactor  will  be  con¬ 
sidered  subsequently o  Let  cp^  be  the  assumed  fast  flux. 


H  ft  M  N  H 


FIGURE  10 


t  i  FOR  A  SLAB  REACTOR  WITH 
*-<V2)  INFINITE  WATER 

REFLECTORS 


D, 

d2 
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.120 

.120 

.275 
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.0115 

.0150 

.0418 

.0926 

qf'  ql2 


vanishing  on  the  boundary  and  normalized  so  that  max  Cpf  “  1°  source  terms 

are  qg  =  2fg(pf  in  Rg  and  =  2fl<pf  in  \  (  2fx  and  2f2  are  the  fast,  ^ 
and  Z2  the  slow  neutron  absorption  cross  sections).  Green* s  function  G  is  the 
solution  of  D1V2G  -  2-jG  which  vanishes  on  the  boundary.  The  flux  (p  is  the 
solution  of  the  integral  equation 

D  (p  =  (D  -Dg)  /<P  fjj  +  (ki2"K22)D2  +  1?1  1* fG  + 

Sg  «2  .  “l  2 

Now  suppose  the  region  extends  to  infinity  and  that  the  source 
terms  are  the  constants  and  Xf£o  Denote  Greenes  function  for  this  problem 

by  G^o  The  flux  distribution,  \}f,  is  the  solution  of 

h*  =  'W  f*  fr  +  <*i2- VH  1°**  +  Sfis It +  ( Zf2‘ 2fl>  {aK 


The  superposed  flux  is  <p^s)  =  i|Kpfs  the  error  to  be  estimated  is  8<p  -  <p  - 
Put  G*  =  G  +  8 Go  Then  we  can  derive  the  following  integral  equation  for  8(p. 

(II. 7)  D  8  cp  =  (DrDg)  /Sep  ~  /gS«P  + 


(II. 8) 


8  Q  =  (Dx-D, z)f Jvjg  -  cpf  ft  g/+  (^i2-^2)  -  %  $*] 

*  Ifi -  *1 1 J*  ^  ^  /°7 

-  », A»rV  /*  ^ir +  M +  |G  +  /s7 
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FINITE  SLAB  REACTOR  WITH  INTERNAL  WATER  REFLECTOR 


As  was  to  be  expected,  §cp  satisfies  an  integral  equation  which 
differs  from  the  equation  satisfied  by  cp  only  in  the  inhomogeneous  term  8  Q. 
Equation  (ll„8)  furnishes  justification  for  superposition  by  multiplication 
since  8  Q  will  be  small  provided  R]_  is  large  enough  for  8  G  to  be  small,  and 
provided  <pf  does  not  vary  rapidly  over  the  regions  which  are  significant  in 
the  integrations .  Hence  R^  should  not  be  too  large  nor  too  near  the  boundary 

of  R^o 

By  way  of  illustrating  what  can  be  expected  of  multiplicative  super- 
position,  the  flux  pattern  created  by  a  slab  water  hole  in  a  finite  slab  reactor 
having  the  STR  constants  lias  been  examined  exactly  and  by  superposition »  (See 
sketch  on  page  54 o) 

We  are  concerned  here  with  a  comparatively  wide  water  hole  very  close 
to  the  boundary  of  the  reactor.  Using  the  STR  constants  given  on  page  25,  and 
putting  <pf  =  cos  ~  the  thermal  flux  is  found  to  be 


•mr  -8  +„4471x  ,  -„4471x 

=  18.2  cos  H  -  6.32x10  e  +  .103  e  , 

=  102.9  009  f  -  120  e-5471*  -  .0218  e+'M71X  , 

=  18.2  oos  H  +  2766  e— 4471*  -  1.692xl0~3  e-4471*. 


-16  <  x  <  10, 


10  <  x  <  14, 


14  <  x  <  16 


The  superposition  flux  is 


=  |-53o3  cosh^"~  o  347  (x— 12)7  102  o  9^  5*?  > 

-  {  44 08  e-°447(x“12)  +  1802}  COS  —  , 


10  <  x  <  14, 


-16  <  x  <  10 
14  <  x  <  16 
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PLICATIVE  SUPERPOSITION  FOR 
■  SLAB  WATER  HOLE  IN  A 
FINITE  SLAB  REACTOR 


J  9  A  9  “]  xn|J  U0J|n9fy/ 


Distance  In  Centimeters 


The  flux  patterns  are  illustrated  in  Figure  11,  page  56-  Evidently 
the  method  of  superposition  provides  a  reasonable  approximation  to  the  correct 
flux. 


3 .  Combined  Superposition 

Now  let  us  look  at  a  three  region  problem.  Region  is  the  reactor 
core,  Rg  is,  for  example,  a  water  hole  inside  the  core,  and  Rg  is  a  reflector 
which  we  assume  extends  to  infinity. 

*  I 


A  method  of  approximating  the  flux  distribution,  (p,  consists  of 
multiplying  that  flux  distribution,  ijr,  corresponding  to  a  constant  fast  flux, 
by  an  appropriate  weighting  factor,  <pf .  The  flux  distribution,  \jr,  corresponding 
to  a  constant  fast  flux,  can  in  turn  be  approximated  by  superposing  the  flux 
patterns  caused  by  R^  and  Rg  separately. 

The  flux  cp  is  the  solution  of  (il.l)  with  qp  =  ^f  Xpf •  The  flux  \|r 
satisfies  the  same  integral  equation  with  q^  =  £  f\ .  Then  3  <p  =  (p  -  ijfcp^  is 
the  solution  of 

D±S  cp  =  L(S  <p)  +  L(ijrtpf)  -  (pfLtt)  +  -  <pf  J  I fG  J 

**1+2+3  **1+2+3 

and  8  (p  will  be  small  if  cpf  does  not  vary  too  rapidly. 
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4o  Conclusion 


In  retrospect  it  can  be  said  that  empirical  evidence  justifies  the 
method  of  superposition  under  circumstances  which  at  first  sight  seem  drastic „ 

In  addition,  there  exist  analytical  expressions  which  provide  a  theoretical 
basis  for  superposition  On  the  other  hand,  attempts  to  estimate  the  super¬ 
position  error  by  the  formulae  developed  here  would  usually  not  be  worth  the 
efforto 

If  it  is  considered  necessary  to  calculate  the  flux  distribution  in  a 
complicated  configuration  more  accurately  than  can  be  done  by  superposition,  then 
it  appears  advantageous  to  solve  the  integral  equation  not  for  the  flux  (p,  but 
for  the  error  8<p  since  this  error  is  generally  small  and  a  slowly  varying  func¬ 
tion  As  was  pointed  out  above,  to  a  first  approximation  D^S  (p  =  SQ,  a  known 
function o 
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IIIo  OTHER  APPLICATIONS 


In  section  1  of  part  III  we  shall  discuss  the  application  of 
integral  equations  to  control  rods,  and  examine  the  superposition  principle 
for  control  rods  and  for  a  control  rod  and  a  diffusing  region0  The  extension 
of  the  integral  equation  method  to  raultigroup  and  characteristic  value  problems 
will  be  outlined  in  section  2o  Other  aspects  of  the  integral  equations — per¬ 
turbation  formulae,  variational  principles,  etc., — will  not*  be  discussed  in 
this  reporto  These  matters  are  not  difficult  although  the  presence  of  surface 
integrals  differentiates  the  integral  equations  from  those  commonly  treated 
in  textbooks o 

1 o  Control  Rods 

The  integral  equation  defining  the  flux  in  the  neighborhood  of  a 
control  rod  is  easily  derived 0  Region  R^  in  the  sketch  on  page  59  is  the  con¬ 
trol  rod,  R^  is  the  diffusing  mediunu  From  the  equations 

p 

DV  (p  -  Sep  +  q  =  0 
dV2G  -  lG  =0 

is  derived  the  integral  equation 


A  The  integral  equation  (ill.l)  has  been  previously  derived  by  M„  Danzker  of 
the  Physics  Department,  Westinghouse  Atomic  Power  Division. 
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The  usual  boundary  condition  is  a  relation  between  <p  and  its  normal  derivative 

_  _  <e 

9n  X  ’ 


where  Us  a  function  of  position  on  the  surface  S-^.  With  this  boundary  con¬ 
dition  the  equation  takes  its  final  form: 


(III.1) 


=  -  S' 


M*M  ■ 

\  Rl+2 


A  detailed  discussion  of  equation  (III.1)  is  unnecessary.  ¥e  point 
out  only  that  the  integral  equation  can  first  be  solved  for  the  flux  on  the  sur¬ 
face  and  then  the  flux  elsewhere  is  obtained  by  integration. 

Zo  Superposition  of  Control  Rods 

The  integral  equation  satisfied  by  the  flux  distribution  around  two 
control  rods  in  an  infinite  reactor  (see  sketch  on  page  61)  can  be  written 


=  ,<p-  /GC  +  JaD 

2+3  *2+3  *1+2+3 


(2)  (3) 

The  flux  patterns  <p  '  and  cp  ,  created  by  the  rods  R^  and  separately, 
satisfy  the  equations 


(i)  -  -  /<*?  *  /< 


^+2+3 


i  =  1,2 


The  superposed  flux  is  defined  to  be  <p^  =  <p^  +  <p^  -  q/£ 


Hence  the  error  S  <p  =  <p  —  '  is  the  solution  of 
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Sep  = 


5G  G  \  c\ 

to+5i  >8* 


/<f +!  >/V5>-  fJ  -  /< 


f ♦ f  fj. 


As  was  the  case  when  R^  and  Rg  are  diffusing  media,  the  error  in  superposition 
depends  on  the  deviations  of  the  individual  flux  patterns  from  their  asymptotic 
values. 

The  last  example  we  shall  mention  is  the  superposition  of  a  control 
rod,  Rg,  and  a  diffusing  region,  Rg,  in  an  infinite  reactor,  R^.  The  integral 
equations  for  (p  and  the  superposition  error  9  are: 

V  -  +  k/)D£  fa 

S2  *2 


/(  to  +  X  )(p  +  /qlG  +  /q2G  ’ 
S  Rq  r2 


i  =  1,2 


DiS  9  =  (VV  +  (k12"  K22)D2  /G8(p 

S2  R2 

-  /(to+\)S<P+  (W  /^"(p(3)~  2^  3  to 


♦  SS  //V5)-  &  -7  f  -  /( I  ♦  S  »^<2)-  £-7 


3o  Generalizations 

The  extension  of  the  integral  equations  derived  above  to  include 
several  energy  groups  of  neutrons  or  to  characteristic  value  problems  offers  no 
difficulty o  Corresponding  to  each  energy  group  there  is  an  integral  equation 
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whose  source  term  is  the  source  of  neutrons  appropriate  to  the  group„  For 
each  group  there  will  exist  a  Green's  function  which  depends  upon  k  and  the 
boundary  conditionso  Apart  from  questions  of  convergence ,  the  characteristic 
value  problem  can  be  attacked  by  iteration  since  the  ratios  of  successive 
iterates  tend  to  a  limit  which  determines  the  reactivity  of  the  pile„  For 
complicated  geometries  iteration  of  the  integral  equation  or  application  of  a 
variational  principle  provides  perhaps  the  most  useful  approach,,  Of  course  the 
variational  technique  can  be  applied  not  only  to  the  integral  equations  but  to 
the  multi-group  differential  equations  as  well,*  A  disadvantage  of  the  varia¬ 
tional  method  is  that  not  only  do  the  fluxes  have  to  be  approximated  but  the 
adjoint  fluxes  also  enter  into  the  equations , 
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